#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <string.h>
#include <inttypes.h>
#include <regex.h>
#include <math.h>
#include <time.h>


// Baffon's Needle
static double d = 4; // distance of parallel lines
static double l = 3; // length of the needle
static int n = 10;    // 8-12, divide first quadrant angular equally into 2^n
static double lambda = 1.96; // alpha = 0.05
static double probabi;
static double mcpi;
static double sigma;
static double error;

// 64 bit ffff
uint64_t x_ffff = 0xffffffffffffffff;

// RNG random number generator
static int p;
static uint64_t s[16];

static inline uint64_t rotl(const uint64_t x, int k) {
	return (x << k) | (x >> (64 - k));}

uint64_t next(void) {
	const int q = p;
	const uint64_t s0 = s[p = (p + 1) & 15];
	uint64_t s15 = s[q];
	const uint64_t result = rotl(s0 * 5, 7) * 9;
	s15 ^= s0;
	s[q] = rotl(s0, 25) ^ s15 ^ (s15 << 27);
	s[p] = rotl(s15, 36);
	return result;}
// end RNG


// rng_init
static int rng_init(){
    printf("Initializing RNG...\n");
    srand((unsigned)time(NULL)); // seeding rand()
    
    p = (rand()%491)+1; // seeding p, 1<= p <= 491 
    printf("Seeding p with %d\n",p);
    
    for (int i = 0; i < 16; i++) {
        s[i] = rand();
        printf("Seeding s[%d] with 0x%x\n",i,s[i]);} // seeding s[]
    
    int mix = 65537; // 65537,491
    for (int i = 1; i <= mix; i++) {
        next(); // mixing
        printf("Initializing XorShift...[%d/%d]\r",i,mix);
        fflush(stdout);}

    printf("\nInitialized RNG successfully.\n");
    return 0;}
// end rng_init


// pow_int
static int pow_int(int a, int b){
    int c = 1;
    for (int i=0; i<b; i++){
        c *= a;}
    return c;}
// end pow_int


// rad_list
static int n2;
static int n2_2;
static double s_list[65536];
static double v_list[65536][2];

static int rad_list(){
    int i; int j; 
    int j_step;
    double m_v;
    
    n2 = pow_int(2,n);
    n2_2 = n2 * 2;
    
    for (i = 0; i < 65536; i++){
        v_list[i][0] = 0; v_list[i][1] = 0;
        s_list[i] = 0;}
    v_list[0][0]=1; v_list[0][1]=0; // x, (1,0)
    v_list[n2][0]=0; v_list[n2][1]=1; // y, (0,1)

    j_step = n2; 
    m_v = 1;
    for (i = 1; i <= n; i++){
        j_step /= 2; 
        j = j_step;
        for ( ; j < n2; j += j_step){
            if (v_list[j][0]==0 && v_list[j][1]==0){

                v_list[j][0] = v_list[j+j_step][0] + v_list[j-j_step][0];
                v_list[j][1] = v_list[j+j_step][1] + v_list[j-j_step][1];
                
                m_v = sqrt( pow(v_list[j][0],2) + pow(v_list[j][1],2) );

                v_list[j][0] /= m_v;
                v_list[j][1] /= m_v;
                }}}

    for (i = 0; i < n2; i++){
        s_list[i] = v_list[i][1];}
    for (i = n2; i< n2_2; i++){
        s_list[i] = v_list[i-n2][0];}
    return 0;}
// end rad_list


// drop_count
uint64_t drop_num = 65536; // by default
uint64_t count_num = 0;

static uint64_t drop_count(uint64_t drop_n){
    uint64_t count_n = 0;
    double x_one = 0;
    int x_rad = 0;
    double d_rand, s_rand;
    for (uint64_t i = 0; i < drop_n; i++){
        x_one = (double) next() / x_ffff;
        x_rad = next() % n2_2;
        d_rand = d * x_one;
        s_rand = l * s_list[x_rad];
        if (s_rand + d_rand > d){
            count_n++;}
        printf("Dropping Needles...[%llu/%llu]\r",i+1,drop_n);
        fflush(stdout);}
    printf("\nDropping Needles Successfully Finished.\n");
    return count_n;}
// end drop_count


// main
int main(int argc, char *argv[]) {
    printf("Baffon's Needle Experiment by Z.X.Sun.\n");

    if (argc == 2) {
        drop_num = (uint64_t)strtoull(argv[1], NULL, 0);
        printf("Drop the Needle %llu times (by user)\n",drop_num);}
    else if (argc == 1){
        printf("Drop the Needle %llu times (by default)\n",drop_num);}
    else {
        printf("Usage:    ");
        printf("buffon.run [number]\n");}

    rng_init();

    rad_list(); 

    count_num = drop_count(drop_num);
    probabi = (double)count_num / drop_num;
    mcpi = (2*l)/(probabi*d);
    sigma = sqrt( (double)count_num / drop_num
                - pow(((double)count_num / drop_num),2) );
    error = lambda * sigma 
                / sqrt((double)drop_num);

    printf("Needle Length: %f\n",l);
    printf("Lines Distance: %f\n",d);
    printf("Drops: %llu\n",drop_num);
    printf("Count: %llu\n",count_num);
    printf("Pi: %.16f\n",mcpi);
    printf("sigma: %.16f\n",sigma);
    printf("error: %.16f\n",error);
    return 0;}
// end main
